1 Introduction

The aim of this notebook is to compare the annotations between slan+Mono/macro/DC vs slan+ cells.

2 Pre-processing

2.1 Load packages

library(Seurat)
library(tidyverse)
library(harmony)
library(here)
library(glue)
library(pheatmap2)

2.2 Parameters

# Colors
color_palette <-  c("#1CFFCE", "#90AD1C", "#C075A6", "#85660D", "#5A5156", "#AA0DFE",   
                    "#F8A19F", "#F7E1A0", "#1C8356", "#FEAF16", "#822E1C", "#C4451C",   
                    "#1CBE4F", "#325A9B", "#F6222E", "#FE00FA", "#FBE426", "#16FF32", 
                    "black",   "#3283FE", "#B00068", "#DEA0FD", "#B10DA1", "#E4E1E3",   
                    "#90AD1C", "#FE00FA", "#85660D", "#3B00FB", "#822E1C", "coral2", 
                    "#1CFFCE", "#1CBE4F", "#3283FE", "#FBE426", "#F7E1A0", "#325A9B",   
                    "#2ED9FF", "#B5EFB5", "#5A5156", "#DEA0FD", "#FEAF16", "#683B79",   
                    "#B10DA1", "#1C7F93", "#F8A19F", "dark orange", "#FEAF16", "#FBE426",  
                    "Brown")

# Functions
process_seurat <- function(seurat_obj, group_by_vars = "gem_id") {
  seurat_obj <- seurat_obj %>%
    FindVariableFeatures() %>%
    ScaleData() %>%
    RunPCA() %>%
    RunHarmony(group.by.vars = group_by_vars) %>%
    RunUMAP(reduction = "harmony", dims = 1:20)
  seurat_obj
}

2.3 Read data

slan_cells <- readRDS(here("scRNA-seq/results/R_objects/slancytes_revision/3.1.seurat_slan_cells_tonsil_classified_deep_annotation_experiments_2_and_3.rds"))
unfiltered <- readRDS(here("scRNA-seq/results/R_objects/slancytes_revision/1.3.seurat_slan_cells_tonsil_filtered_experiments_2_and_3.rds"))
table(slan_cells$gem_id)
## 
## q0qp6qyu_nkrug5p0 rssxkrxg_edqhwe1o 
##               914               741

3 Dimensionality reduction

slan_cells <- process_seurat(slan_cells)
slan_cells$library_name2 <- str_remove(slan_cells$library_name, "BCLL_10_")
DimPlot(slan_cells, group.by = "library_name2") +
  theme(
    legend.text = element_text(size = 15),
    legend.position = "top",
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    plot.title = element_blank()
  )

slan_markers <- read.delim(here("data/raw_data_figures/BA2020Faseb_SLANmarkers.txt"))
sel_genes <- slan_markers$Gene[slan_markers$Gene %in% rownames(slan_cells)]
slan_cells <- AddModuleScore(
  slan_cells,
  features = list(
    SLAN = slan_markers$Gene[slan_markers$CellType == "SLAN"],
    DC = slan_markers$Gene[slan_markers$CellType == "DC"],
    MAC = slan_markers$Gene[slan_markers$CellType == "MAC"]
  ),
  name = c("SLAN", "DC", "MAC")
)
FeaturePlot(slan_cells, c("SLAN1", "DC2", "MAC3")) &
  scale_color_viridis_c(option = "magma") &
    theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
  )

DimPlot(slan_cells, group.by = "annotation_20220619", cols = color_palette)

FeaturePlot(slan_cells, "annotation_probability") &
  scale_color_viridis_c(option = "magma")

DimPlot(slan_cells, group.by = "annotation_20220619", cols = color_palette)

slan_cells@meta.data %>%
  group_by(annotation_20220619, library_name) %>%
  summarize(n_cells = n()) %>%
  group_by(library_name) %>%
  mutate(total_cells = sum(n_cells), pct_cells = n_cells / total_cells * 100) %>%
  ggplot(aes(library_name, pct_cells, fill = annotation_20220619)) +
    geom_col() +
    labs(y = "Percentage of cells (%)") +
    scale_fill_manual(values = color_palette) +
    theme_classic() +
    theme(axis.title.x = element_blank(),
          legend.title = element_blank())

VlnPlot(slan_cells, c("SLAN1", "DC2", "MAC3"), group.by = "library_name")

Remove lingering poor-quality cells/doublets

slan_cells$doublet_score <- unfiltered$doublet_score[colnames(slan_cells)]
FeaturePlot(slan_cells, features = "doublet_score") &
  scale_color_viridis_c(option = "magma")

slan_cells <- FindNeighbors(slan_cells, reduction = "harmony", dims = 1:20, k.param = 15)
slan_cells <- FindClusters(slan_cells, resolution = 1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1655
## Number of edges: 51500
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7497
## Number of communities: 14
## Elapsed time: 0 seconds
markers_9 <- FindMarkers(slan_cells, ident.1 = "9", ident.2 = "1", only.pos = TRUE, logfc.threshold = 0.75)
DT::datatable(markers_9, options = list(scrollX = TRUE))
markers_3 <- FindMarkers(slan_cells, ident.1 = "3", only.pos = TRUE, logfc.threshold = 0.75)
DT::datatable(markers_3, options = list(scrollX = TRUE))
DimPlot(slan_cells)

VlnPlot(slan_cells, "nFeature_RNA", pt.size = 0) & NoLegend()

slan_cells <- slan_cells[, slan_cells$seurat_clusters != "9"]
DimPlot(slan_cells)

slan_cells <- process_seurat(slan_cells)
DimPlot(slan_cells)

4 Subcluster

slan_cells <- FindNeighbors(slan_cells, reduction = "harmony", dims = 1:20)
slan_cells <- FindClusters(slan_cells, resolution = 0.15)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1571
## Number of edges: 54794
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9138
## Number of communities: 5
## Elapsed time: 0 seconds
DimPlot(slan_cells)

slan_sub <- slan_cells[, slan_cells$seurat_clusters == "1" & slan_cells$library_name2 == "SLANCYTES"]
slan_sub <- slan_sub %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA() %>%
  # RunHarmony(group.by.vars = "gem_id") %>%
  RunUMAP(reduction = "pca", dims = 1:20)
DimPlot(slan_sub, group.by = "gem_id")

FeaturePlot(slan_sub, c("SELENOP", "C1QA", "MMP9", "MAF"))

# Heatmap
goi_slan <- c("SELENOP", "APOC1", "APOE", "FCGR2A", "MMP9", "MMP12",
              "C1QA", "C1QB", "HLA-DRA", "HLA-DRB1")
slan_sub <- FindNeighbors(slan_sub, reduction = "pca", dims = 1:20, k.param = 7)
slan_sub <- FindClusters(slan_sub, resolution = 1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 370
## Number of edges: 8728
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.5370
## Number of communities: 6
## Elapsed time: 0 seconds
DimPlot(slan_sub)

slan_sub$subcluster <- case_when(
  slan_sub$seurat_clusters %in% c("0", "1") ~ "0",
  slan_sub$seurat_clusters %in% c("3", "4") ~ "1",
  slan_sub$seurat_clusters == "2" ~ "2",
  slan_sub$seurat_clusters == "5" ~ "3"
)
barcodes <- slan_sub@meta.data %>%
  rownames_to_column("barcode") %>%
  arrange(subcluster) %>%
  pull("barcode")
input_mat <- as.matrix(slan_sub[["RNA"]]@data[goi_slan, barcodes])
pheatmap2(
  input_mat,
  scale = "none",
  show_rownames = TRUE,
  show_colnames = FALSE, 
  border_color = NA,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  fontsize_row = 6
)

5 Save

saveRDS(slan_cells, here("scRNA-seq/results/R_objects/slancytes_revision/4.1.seurat_slan_cells_tonsil_dimensionality_reduction_experiments_2_and_3.rds"))
saveRDS(slan_sub, here("scRNA-seq/results/R_objects/slancytes_revision/4.2.seurat_slan_cells_tonsil_subset_slancytes_experiments_2_and_3.rds"))

6 Session Information

sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pheatmap2_0.1.0    fastcluster_1.2.3  parallelDist_0.2.6 gtable_0.3.1       scales_1.2.1       RColorBrewer_1.1-3 glue_1.6.2         here_1.0.1         harmony_0.1.1      Rcpp_1.0.10        forcats_1.0.0      stringr_1.5.0      dplyr_1.1.0        purrr_1.0.1        readr_2.1.3        tidyr_1.3.0        tibble_3.1.8       ggplot2_3.4.0      tidyverse_1.3.2    SeuratObject_4.1.3 Seurat_4.3.0       BiocStyle_2.26.0  
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.1           backports_1.4.1        plyr_1.8.8             igraph_1.3.5           lazyeval_0.2.2         sp_1.6-0               splines_4.2.0          crosstalk_1.2.0        listenv_0.9.0          scattermore_0.8        digest_0.6.31          htmltools_0.5.4        fansi_1.0.4            magrittr_2.0.3         tensor_1.5             googlesheets4_1.0.1    cluster_2.1.4          ROCR_1.0-11            limma_3.54.1           tzdb_0.3.0             globals_0.16.2         modelr_0.1.10          RcppParallel_5.1.6     matrixStats_0.63.0     timechange_0.2.0       spatstat.sparse_3.0-0  colorspace_2.1-0       rvest_1.0.3            ggrepel_0.9.3          haven_2.5.1            xfun_0.37              crayon_1.5.2           jsonlite_1.8.4         progressr_0.13.0       spatstat.data_3.0-0    survival_3.5-0         zoo_1.8-11             polyclip_1.10-4        gargle_1.3.0           leiden_0.4.3           future.apply_1.10.0    abind_1.4-5            DBI_1.1.3              spatstat.random_3.1-3  miniUI_0.1.1.1         viridisLite_0.4.1      xtable_1.8-4           reticulate_1.28        DT_0.27                htmlwidgets_1.6.1      httr_1.4.4            
##  [52] ellipsis_0.3.2         ica_1.0-3              farver_2.1.1           pkgconfig_2.0.3        sass_0.4.5             uwot_0.1.14            dbplyr_2.3.0           deldir_1.0-6           utf8_1.2.3             labeling_0.4.2         tidyselect_1.2.0       rlang_1.0.6            reshape2_1.4.4         later_1.3.0            munsell_0.5.0          cellranger_1.1.0       tools_4.2.0            cachem_1.0.6           cli_3.6.0              generics_0.1.3         broom_1.0.3            ggridges_0.5.4         evaluate_0.20          fastmap_1.1.0          yaml_2.3.7             goftest_1.2-3          knitr_1.42             fs_1.6.0               fitdistrplus_1.1-8     RANN_2.6.1             pbapply_1.7-0          future_1.31.0          nlme_3.1-162           mime_0.12              ggrastr_1.0.1          xml2_1.3.3             compiler_4.2.0         rstudioapi_0.14        beeswarm_0.4.0         plotly_4.10.1          png_0.1-8              spatstat.utils_3.0-1   reprex_2.0.2           bslib_0.4.2            stringi_1.7.12         highr_0.10             lattice_0.20-45        Matrix_1.5-3           vctrs_0.5.2            pillar_1.8.1           lifecycle_1.0.3       
## [103] BiocManager_1.30.19    spatstat.geom_3.0-6    lmtest_0.9-40          jquerylib_0.1.4        RcppAnnoy_0.0.20       data.table_1.14.6      cowplot_1.1.1          irlba_2.3.5.1          httpuv_1.6.8           patchwork_1.1.2        R6_2.5.1               bookdown_0.32          promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3          vipor_0.4.5            parallelly_1.34.0      codetools_0.2-19       MASS_7.3-58.2          assertthat_0.2.1       rprojroot_2.0.3        withr_2.5.0            sctransform_0.3.5      parallel_4.2.0         hms_1.1.2              rmarkdown_2.20         googledrive_2.0.0      Rtsne_0.16             spatstat.explore_3.0-6 shiny_1.7.4            lubridate_1.9.1        ggbeeswarm_0.7.1